#define _USE_MATH_DEFINES

#include <iostream>
#include <fstream>
#include <datatable.h>
#include <bspline.h>
#include <bsplinebuilder.h>
#include <nlohmann/json.hpp>
#include <vector>
#include <cmath>

#include "polynomial.h"

using std::cout;
using std::endl;

using namespace SPLINTER;

// Six-hump camelback function
double f(DenseVector x)
{
    assert(x.rows() == 2);
    return (4 - 2.1 * x(0) * x(0) + (1 / 3.) * x(0) * x(0) * x(0) * x(0)) * x(0) * x(0) + x(0) * x(1) + (-4 + 4 * x(1) * x(1)) * x(1) * x(1);
}

int main(int argc, char *argv[])
{
    std::cout << 123 << std::endl;
    Polynomial aa(4, 1.0, 2.0, 3.0, 4.0);
    Polynomial bb(2, 2.0, 1.0);
    DenseVector a = aa.getCoefficients(), b = bb.getCoefficients();

    DenseVector d(a);
    std::cout << "aa:" << aa.to_string() << std::endl;
    std::cout << "bb:" << bb.to_string() << std::endl;
    Polynomial cc = aa.derivation().derivation().derivation().derivation();
    std::cout << "123" << std::endl;
    std::cout << "cc:" << cc.to_string() << std::endl;

    for (int i = 0; i < a.size(); i++)
    {
        std::cout << a[i] << ",";
    }
    std::cout << std::endl;
    for (int i = 0; i < b.size(); i++)
    {
        std::cout << b[i] << ",";
    }
    std::cout << std::endl;
    for (int i = 0; i < d.size(); i++)
    {
        std::cout << d[i] << ",";
    }
    std::cout << std::endl;
    nlohmann::json j;
    {
        std::ifstream file("input.json");
        j = nlohmann::json::parse(file);
        file.close();
    }

    DataTable samples0;
    { // Create new DataTable to manage samples
        DataTable samples;

        // Sample the function
        DenseVector x(2);
        double y;
        for (int i = 0; i < 5; i++)
        {
            for (int j = 0; j < 5; j++)
            {
                // Sample function at x
                x(0) = i * 0.1;
                x(1) = j * 0.1;
                y = f(x);

                // Store sample
                samples.addSample(x, y);
                // std::cout<<"["<<x(0)<<","<<x(1)<<","<<y<<"],"<<std::endl;
            }
        }
        samples0 = samples;
    }

    auto data = j["data"];

    for (auto it = data.begin(); it != data.end(); it++)
    {
        std::vector<std::string> names = (*it)["names"];
        auto samples = (*it)["samples"];
        auto targets = (*it)["targets"];
        DataTable samples_table;
        auto it_samples = samples.begin();
        size_t y_index = 0;
        for (int i = 0; i < names.size(); i++)
        {
            if (!strcmp(names[i].c_str(), "y"))
            {
                y_index = i;
                break;
            }
        }
        DenseVector x(names.size() - 1);
        double y;
        if (it_samples->size() == names.size())
        {

            for (; it_samples != samples.end(); it_samples++)
            {
                int x_index = 0;
                for (int i = 0; i < names.size(); i++)
                {
                    if (i != y_index)
                        x(x_index++) = (double)(*it_samples)[i];
                    else
                        y = (*it_samples)[i];
                }
                samples_table.addSample(x, y);
            }
        }
        else if (samples.size() == names.size())
        {
            int len = samples[0].size();
            for (int i = 1; i < samples.size(); i++)
            {
                if (samples[i].size() != len)
                {
                    throw Exception("main: Input type error.");
                }
            }
            for (int j = 0; j < len; j++)
            {
                int x_index = 0;
                for (int i = 0; i < names.size(); i++)
                {
                    if (i != y_index)
                        x(x_index++) = samples[i][j];
                    else
                        y = samples[i][j];
                }
                samples_table.addSample(x, y);
            }
        }
        else
        {
            throw Exception("main: Input type error.");
        }

        // Build B-splines that interpolate the samples
        BSpline bspline1 = BSpline::Builder(samples_table).degree(0).build();
        BSpline bspline3 = BSpline::Builder(samples_table).degree(3).build();

        // Build penalized B-spline (P-spline) that smooths the samples
        BSpline pspline = BSpline::Builder(samples_table)
                              .degree(3)
                              .smoothing(BSpline::Smoothing::PSPLINE)
                              .alpha(0.03)
                              .build();

        /* Evaluate the approximants at x = (1,1)
         * Note that the error will be 0 at that point (except for the P-spline, which may introduce an error
         * in favor of a smooth approximation) because it is a point we sampled at.
         */
        auto it_targets = targets.begin();
        if (it_targets->size() == names.size() - 1)
        {
            for (; it_targets != targets.end(); it_targets++)
            {
                for (int i = 0; i < it_targets->size(); i++)
                {
                    x(i) = (*it_targets)[i];
                }
                cout << "-----------------------------------------------------" << endl;
                cout << "Function at x:                 " << f(x) << endl;
                cout << "Linear B-spline at x:          " << bspline1.eval(x) << endl;
                cout << "Cubic B-spline at x:           " << bspline3.eval(x) << endl;
                cout << "P-spline at x:                 " << pspline.eval(x) << endl;
                cout << "-----------------------------------------------------" << endl;
            }
        }
    }

    return 0;
}